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ABSTRACT 

Estimates of the Galactic coalescence rate (TZ) of close binaries with two 
neutron stars (NS-NS) are known to be uncertain by large factors (about two 
orders of magnitude) mainly due to the small number of systems detected as 
binary radio pulsars. We present an analysis method that allows us to estimate 
the Galactic NS-NS coalescence rate using the current observed sample and, 
importantly, to assign a statistical significance to these estimates and to calculate 
the allowed ranges of values at various confidence levels. The method involves 
the simulation of selection effects inherent in all relevant radio pulsar surveys and 
a Bayesian statistical analysis for the probability distribution of the rate. The 
most likely values for the total Galactic coalescence rate (7£ pea k) lie in the range 
2 — 60 Myr -1 depending on different pulsar population models. For our reference 
model 1, where the most likely estimates of pulsar population properties are 
adopted, we obtain lZ tot = 8+5 Myr -1 at a 68% statistical confidence level. The 
corresponding range of expected detection rates of NS-NS inspiral are 3+2 X 10~ 3 
y r -i f or initial LIGO and 18+^ yr" 1 for the advanced LIGO. 

Subject headings: binaries: close-gravitational waves-stars: neutron 



1. INTRODUCTION 

The detection of the double neutron star (NS-NS) prototype PSR B1913+16 as a binary 
pulsar (Hulse & Taylor 1975) and its orbital decay due to emission of gravitational waves 
(Taylor et al. 1979) has inspired a number of quantitative estimates of the coalescence rate of 
NS-NS binaries (Clark et al. 1979; Narayan et al. 1991; Phinney 1991; Curran & Lorimer 1995). 
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In general, the coalescence rate of NS-NS binaries can be calculated based on: (a) our the- 
oretical understanding of their formation (see Belczynski & Kalogera 2001 for a review and 
application of this approach); (b) the observational properties of the pulsars in the binary 
systems and the modeling of pulsar survey selection effects (see e.g. Narayan 1987). Interest 
in these coalescence derives from an intrinsic motivation of understanding their origin and 
evolution and their connections to other NS binaries. However, significant interest derives 
from their importance as gravitational-wave sources for the upcoming ground-based laser in- 
terferometers (such as LIGO) and their possible association with 7-ray burst events (Popham 
et al. 1998 and references therein). 

The traditional way of calculating the coalescence rate based on observations involves 
an estimate of the scale factor, an indicator for the number of pulsars in our Galaxy 
with the same spin period and luminosity (Narayan 1987). Corrections must then be ap- 
plied to these scale factors to account for the faint end of the pulsar luminosity func- 
tion, the beamed nature of pulsar emission, and uncertainties in the assumed spatial dis- 
tribution. The estimated total number in the Galaxy can then be combined with es- 
timates of their lifetimes to obtain a coalescence rate, 1Z. This method was first ap- 
plied by Narayan et al. (1991) and Phinney (1991) and other investigators who followed 
(Curran & Lorimer 1995; van den Heuvel & Lorimer 1996). Various correction factors were 
(or were not) included at various levels of completeness. Summaries of these earlier studies 
can be found in Arzoumanian et al. (1999) and Kalogera et al. (2001; hereafter KNST). 
The latter authors examined all possible uncertainties in the estimates of the coalescence 
rate of NS-NS binaries in detail, and pointed a small-number bias that introduces a large 
uncertainty (more than two orders of magnitude) in the correction factor for the faint-pulsar 
population that must be applied to the rate estimate. They obtained a total NS-NS rate 
estimate in the range 1Z = 10~ 6 — 5x 10 _4 yr _1 , with the uncertainty dominated by the small- 
number bias. Earlier studies, which made different assumptions about the pulsar properties 
(e.g. luminosity and spatial distributions and lifetimes), are roughly consistent with each 
other (given the large uncertainties). Estimated ranges of values until now were not as- 
sociated with statistical significance statements and an "all-inclusive" estimated Galactic 
coalescence rate lies in the range ~ 10~ 7 — 10~ 5 yr _1 . 

The motivation for this paper is to update the scale factor calculations using the most re- 
cent pulsar surveys, and present a statistical analysis that allows the calculation of statistical 
confidence levels associated with rate estimates. We consider the two binaries found in the 
Galactic disk: PSRB1913+16 (Hulse & Taylor 1975) and PSR B1534+12 (Wolszczan 1991). 
Following the arguments made by Phinney (1991) and KNST, we do not include the globu- 
lar cluster system PSR B2127+11C (Prince et al. 1991); this system will be the subject of a 
later paper. Radio-pulsar-survey selection effects are taken into account in the modeling of 
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pulsar population. As described in what follows, the small-number bias and the effect of a 
luminosity function are implicitly included in our analysis, and therefore a separate correc- 
tion factor is not needed. For each population model of pulsars, we derive the probability 
distribution function of the total Galactic coalescence rate weighted by the two observed 
binary systems. In our results we note a number of important correlations between 7^ pe ak 
and model parameters that are useful in generalizing the method. We extrapolate the Galac- 
tic rate to cover the detection volume of LIGO and estimate the detection rates of NS-NS 
inspiral events for the initial and advanced LIGO. 

The plan for the rest of this paper is as follows. In §2, we describe our analysis method 
in a qualitative way. Full details of the various pulsar population models and survey selection 
effects are then given in §3 and §4 respectively. In §5, we derive the probability distribution 
function for the total Galactic coalescence rate and calculate the detection rate of LIGO. In 
§6, we summarize our results and discuss a number of intriguing correlations between various 
physical quantities. Finally, in §7, we discuss the results and compare them with previous 
studies. 



2. BASIC ANALYSIS METHOD 

Our basic method is one of "forward" analysis. By this we mean that we do not attempt 
to invert the observations to obtain the total number of NS-NS binaries in the Galaxy. 
Instead, using Monte Carlo methods, we populate a model galaxy with NS-NS binaries (that 
match the spin properties of PSR B1913+16 and PSR B1534+12) with pre-set properties in 
terms of their spatial distribution and radio pulsar luminosity function. Details about these 
"physical models" are given in §3. 

For a given physical model, we produce synthetic populations of different total numbers 
of objects (N tot ). We then produce a very large number of Monte Carlo realizations of 
such pulsar populations and determine the number of objects (A b s ) that are observable 
by all large-scale pulsar surveys carried out to date by detailed modeling of the detection 
thresholds of these surveys. This analysis utilizes code to take account of observational 
selection effects in a self-consistent manner, developed and described in detail by Lorimer et 
al. (1993; hereafter LBDH) which we summarize in §4. Performing this analysis for many 
different Monte Carlo realizations of the physical model allows us to examine the distribution 
of Aobs- We find, as expected and assumed by other studies, that this distribution closely 
follows Poisson statistics, and we determine the best-fit value of the mean of the Poisson 
distribution A for each population model and value of A to t- 
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The calculations described so far are performed separately for each of PSR B1913+16 
and PSR B1534+12 so that we obtain separate best-fit A values for the Poisson distribu- 
tions. Doing the analysis in this way allows us to calculate the likelihood of observing just 
one example of each pulsar in the real-world sample. Given the Poissonian nature of the dis- 
tributions this likelihood is simply: P(l; A) = Aexp(— A). We then calculate this likelihood 
for a variety of assumed iV to t values for each physical model. 

The probability distribution of the total coalescence rate 7^ ot is derived using the 
Bayesian analysis and the calculated likelihood for each pulsar (described in detail in §5). 
The derivation of this probability distribution allows us to calculate the most probable rate 
as well as determine its ranges of values at various statistical confidence levels. Finally, we 
extrapolate the Galactic rate to the volume expected to be reached by LIGO and calculate 
the detection rate, 1Z det (see §6). 

3. MODELS FOR THE GALACTIC PULSAR POPULATION 

Our model pulsar populations are characterized by a Galactocentric radius (R), vertical 
distance (Z) from the Galactic plane and radio luminosity (L). Assuming that the distribu- 
tions of each of these parameters are independent, the combined probability density function 
(PDF) of the model pulsar population can be written as: 



where iPr(R), ipz(Z) and (f>(L), are the individual PDFs of R, Z and L, respectively. In all 
models considered, we assume azimuthal symmetry about the Galactic center. 

The spatial distribution of pulsars is rather loosely constrained, but we find that it does 
not affect the results significantly for a wide range of models. For the radial and the vertical 
PDFs, we consider Gaussian and exponential forms with different values of the radial Rq 
and the vertical Z scale. In our reference model, we assume a Gaussian PDF for the radial 
component and an exponential PDF for the vertical component. Hence, the combined spatial 
PDF is given by: 



We set R = 4.0 kpc and Z = 1.5 kpc as standard model parameters. Following Narayan et 
al. (1991), these and other values considered reflect the present-day spatial distribution of the 
NS-NS binary population after kinematic evolution in the Galactic gravitational potential. 

Having assigned a position of each pulsar in our model galaxy, for later computational 
convenience, we store the positions as Cartesian x, y, z coordinates, where the Galactic center 



f{R,Z,L) dRdZdL = ij R (R)2nRdR ip z (Z)dZ </>(L)dL, 



(1) 




(2) 
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is defined as (0.0,0.0,0.0) kpc and the position of the Earth is assumed to lie 8.5 kpc from 
the center along the x — y plane, i.e. (8.5,0.0,0.0). From these definitions, the distance d 
to each pulsar from the Earth can be readily calculated, as well as the apparent Galactic 
coordinates I and b. 

For the luminosity PDF, we follow the results of Cordes & Chernoff (1997) and adopt 
a power-law function of the form 

<P(L) = (p-m-^, (3) 

where L > L min and p > 1. The cut-off luminosity, L min , and the exponent p are the model 
parameters. Cordes & Chernoff (1997) found L min = l.l^ois m Jy kpc 2 and p = 2.0 ± 0.2 
at 68% confidence. We set p = 2.0 and L min = 1.0 mJy kpc 2 for our reference model. 
Throughout this paper, luminosities are defined to be at the observing frequency v = 400 
MHz. 

Having defined the position and luminosity of each pulsar in our model Galaxy, the 
final step in defining the model population is to calculate a number of derived parame- 
ters required to characterize the detection of the model pulsars: dispersion measure (DM), 
scatter-broadening time (r) and sky background temperature (T s k y ). To calculate DM and 
r, we use the software developed by Taylor & Cordes (1993) to integrate their model of 
the free-electron column density along the line of sight to each pulsar defined by its model 
Galactic coordinates I and b out to its distance d. Frequency scaling of r to different survey 
frequencies is done assuming a Kolmogorov turbulence spectrum with a spectral index of 
-4.4 1 . Finally, given the model Galactic coordinates of each pulsar, the sky background noise 
temperature at 408 MHz (T s k y ) is taken from the all-sky catalog of Haslam et al. (1981). 
Scaling T sky to other survey frequencies assumes a spectral index of -2.8 (Lawson et al. 1987). 



4. PULSAR SURVEY SELECTION EFFECTS 

Having created a model pulsar population with a given spatial and luminosity distri- 
bution, we are now in a position to determine the fraction of the total population which 
are actually detectable by current large-scale pulsar surveys. To do this, we need to calcu- 
late, for each model pulsar, the effective signal-to-noise ratio it would have in each survey 
and compare this with the corresponding detection threshold. Only those pulsars which are 



1 Although recent studies suggest a variety of spectral indices for r (Lohmcr et al. 2001), the effects of 
scattering turn out to be negligible in this study since the detections of NS-NS binaries are limited by 
luminosity to nearby systems. 
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nominally above the threshold count as being detectable. After performing this process on 
the entire model pulsar population of size N tot , we are left with a sample of N ohs pulsars that 
are nominally detectable by the surveys. Repeating this process many times, we can deter- 
mine the probability distribution of iV bs which we then use to constrain the population and 
coalescence rate of NS-NS binaries. In this section we discuss our modeling of the various 
selection effects which limit pulsar detection. 



The main factors affecting the signal-to-noise ratio (a) of a pulsar search can be sum- 
marized by the following expression 



where S v is the apparent flux density at the survey frequency u, G is the gain of the telescope, 
T is the effective system noise temperature (which includes a contribution T sky from the sky 
background described in the previous section), P is the pulse period, Av is the observing 
bandwidth, t is the integration time and w e is the pulse width. More exact expressions are 
given in the detailed description of the survey selection effects in §2 of LBDH (in particular 
see their eqns. 14-18) which we adopt in this work. In what follows, we describe the salient 
points relevant to this study. 

For each model pulsar, with known 400-MHz luminosity L and distance d, we calculate 
the apparent 400-MHz flux density S400 = L/d 2 . Since not all pulsar surveys are carried out 
at 400 MHz, we need to scale 6400 to take account of the steep radio flux density spectra of 
pulsars. Using a simple power law of the form S„ oc u a , where a is the spectral index, we 
can calculate the flux S u at any frequency v as: 



Following the results of Lorimer et al. (1995), in all simulations, spectral indices were drawn 
from a Gaussian PDF with a mean of -1.6 and standard deviation 0.4. 

The telescope gain, system noise, bandwidth and integration time are well-known pa- 
rameters for any given survey and the detailed models we use take account of these. In 
addition to the surveys considered by LBDH, we also model surveys listed by Curran & 
Lorimer (1995), and more recent surveys at Green Bank (Sayer, Nice & Taylor 1997) and 
Parkes (Lyne et al. 2000; Manchester et al. 2001; Edwards et al. 2001). A complete list of 
the surveys considered, and the references to the relevant publications is given in Table 1. 



4.1. Survey Parameters 




(4) 




(5) 
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Up to this point in the simulations, the model parameters are identical for both PSRs 
B1913+16 and B1534+12. Since we are interested in the individual contributions each 
of these systems make to the total Galactic merger rate of NS-NS binaries similar to these 
systems, we fix the assumed spin periods P and intrinsic pulse widths w to the values of each 
pulsar and perform separate simulations over all physical models considered. The assumed 
pulse widths are 10 ms and 1.5 ms respectively for PSRs B1913+16 and B1534+12. The 
effective pulse width w e required for the signal-to-noise calculation must take into account 
pulse broadening effects due to the interstellar medium and the response of the observing 
system. The various contributions are summarized by the quadrature sum: 

w 2 = w 2 + r 2 + tg amp + + t 2 ADM , (6) 

where r is the scatter-broadening timescale calculated from the Taylor & Cordes (1993) 
model, t samp is the data sampling interval in the observing system, £dm * s the dispersive 
broadening across an individual frequency channel and £adm i s the pulse broadening due to 
dedispersion at a slightly incorrect dispersion measure. All of these factors are accounted 
for in our model described in detail in LBDH. 



4.2. Doppler Smearing 

For binary pulsars, we need to take account of the reduction in signal-to-noise ratio due 
to the Doppler shift in period during an observation. This was not considered in LBDH since 
their analysis was concerned only with isolated pulsars. For observations of NS-NS binaries, 
however, where the orbital periods are of the order of 10 hours or less, the apparent pulse 
period can change significantly during a search observation causing the received power to be 
spread over a number of frequency bins in the Fourier domain. As all the surveys considered 
in this analysis search for periodicities in the amplitude spectrum of the Fourier transform of 
the time series, a signal spread over several bins can result in a loss of signal-to-noise ratio. 
To take account of this effect in our survey simulations, we need to multiply the apparent 
flux density of each model pulsar by a "degradation factor", F. 

To calculate the appropriate F values to use, we generate synthetic pulsar search data 
containing signals with periods and duty cycles similar to PSR B1913+16 and PSR B1534+12. 
These data are then passed through a real pulsar search code which is similar to those in use 
in the large-scale surveys (see Lorimer et al. 2000 for details). For each of the two pulsars, 
we first generate a control time series in which the signal has a constant period and find 
the resulting signal-to- noise ratio, <7 con troi, reported by the search code. We then generate 
a time series in which the pulses have identical intensity but are modulated in period ac- 
cording to the appropriate orbital parameters of each binary pulsar. From the resulting 
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search signal-to-noise ratio, aginary, the degradation factor F = aginary /^control- Significant 
degradation occurs, therefore, when F<1. Since accumulated Doppler shift, and therefore 
F, is a strong function of the orbital phase at the start of a given observation, for both 
binary systems, we calculate the mean value of F for a variety of starting orbital phases 
appropriately weighted by the time spent in that particular part of the orbit. 

A similar analysis was made by Camilo et al. (2000) for the millisecond pulsars in 
47 Tucanae. In this paper, where we are interested in the degradation as a function of 
integration time, we generate time series with a variety of lengths between 1 minute and 1 
hour using sampling intervals similar to those of the actual surveys listed in Table 1. The 
results are summarized in Figure 1, where we plot average F versus integration time for 
both sets of orbital parameters. As expected, surveys with the longest integration times 
are most affected by Doppler smearing. For the Parkes Multibeam survey (Lyne et al. 2000; 
Manchester et al. 2001), which has an integration time of 35 min, mean values of F are 0.7 
and 0.3 for PSR B1913+16 and PSR B1534+12 respectively 2 . The greater degradation for 
PSR B1534+12 is due to its mildly eccentric orbit (e ~ 0.3 versus 0.6 for PSR B1913+16) 
which results in a much more persistent change in apparent pulse period when averaged 
over the entire orbit. For the Jodrell Bank and Swinburne surveys (Nicastro et al. 1995; 
Edwards et al. 2001), which both have integration times of order 5 min, we find F ~ 0.9 for 
both systems. For all other surveys, which have significantly shorter integration times, no 
significant degradation is seen, and we take F — 1. 



5. STATISTICAL ANALYSIS 

In this section we describe in detail the derivation of the probability distribution of 
the Galactic coalescence rate TZ. The analysis method makes use of Bayesian statistics and 
takes into account the rate contributions of both observed NS-NS binaries. At the end of 
the section we derive the associated detection rates for LIGO. 



2 In order to improve on the sensitivity to binary pulsars, the Parkes Multibeam survey data are now 
being reprocessed using various algorithms designed to account for binary motion during the integration 
time (Faulkner et al. 2002) 
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5.1. The Rate Probability Distribution for Each Observed NS-NS Binary 

As already mentioned in § 2, for each of the two observed NS-NS binaries (PSR B1913+16 
or PSR B1534+12) we generate pulsar populations in physical and radio luminosity space 
with pulse periods and widths fixed to the observed ones and with different absolute nor- 
malizations, i.e., total number N tot of pulsars in the Galaxy. We generate large numbers 
of "observed" pulsar samples by modeling the pulsar survey selection effects (see §4) and 
applying them on these model populations of PSR B1913+i6-like and PSR B1534+12-like 
pulsars separately (see §3). For a fixed value of N tot , we use these "observed" samples to 
calculate the distribution of the number of objects in the samples. One might expect that 
the number of observed pulsars in a sample N ohs follows very closely a Poisson distribution: 

\N ohs -X 

P(N ohs ;X) = , (7) 

-< >obs- 

where, by definition, A =< iV bs >• We confirm our expectation by obtaining excellent 
formal fits to the Monte Carlo data using such a distribution and calculating the best-fitting 
value of A. We vary N tot in the range 10 — 10 4 and find that A is linearly correlated with 
N tot : 

A = aN tot , (8) 

where a is a constant that depends on the properties (space and luminosity distributions 
and pulse period and width) of the Galactic pulsar population. Examples of the Poisson 
fits and the linear correlations are shown in Figures 2 and 3, respectively, for our reference 
model 1 (see Table 2 and §6). 

The main step in deriving a rate probability distribution for each of the observed systems 
is to first derive the probability distribution of the total number of pulsars like the observed 
ones in the Galaxy. We obtain the latter by applying Bayes' theorem: 

P(H\DX) = P(H\X) P ^ HX \ (9) 
vi ) P{D\X) ' v ; 

where P(H\DX) is the probability of a model hypothesis H given data D and model priors 
X, P(D\HX) is the likelihood of the data given a model hypothesis and priors, P(H\X) is 
the probability of a model hypothesis in the absence of any data information, and P(D\X) 
is the model prior probability, which acts as a normalization constant. 

In the present work, we make following identifications: 
D : is the real observed sample 
H : is A proportional to N tot 

X : is the population model (space and luminosity distributions and pulse period and width) 
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With these identifications, P(D\HX) is the likelihood of the real observed sample (one 
"PSR B1913+16-like" and one "PSR B1534+12-like" pulsar) and is obtained by the best- 
fitting Poisson distribution: 

P(D\HX) = P(l;X(N tot ),X) = X(N tot )e~ x ^\ (10) 

In the absence of any data information, the absolute normalization of the model population, 
i.e., the total pulsar number A'tot and hence A is expected to be independent of the shape of 
the population distributions and properties represented by X. Therefore, the probability of 
A given a set of assumptions X for the model Galactic population is expected to be flat: 

P{H\X) = P(X(N tot )\X) = constant, (11) 

and is essentially absorbed by the model prior probability P(D\X) as a normalization con- 
stant. The probability distribution of A then is given by: 

P(X\DX) = ^P(1;A,I) = constant x P(1;A,X). (12) 

We impose the normalization constraint: J °° P(X\DX)dX = 1 and find that P(H\X)/P(D\X) 
1 and 

P(X\DX) = P(1;X,X) = Ae~ A . (13) 

Note that based on the above expression, the maximum value of P(X) equal to e _1 always 
occurs at A = 1 or at N tot = a" 1 (see eq.(8)). It is straightforward to calculate P(N tot ): 



P(N tot ) = P(X) 



dX 



dN tot 

a 2 N tot e~ aN ^ (14) 



For a given total number of pulsars in the Galaxy, we can calculate their rate using 
estimates of the associated pulsar beaming correction factor / b and lifetime ri ifc : 

K = -*±f h . (15) 

Tlifc 

Note that this estimate is equivalent to previous studies (KNST and references therein) 
where the concept of the scale factor is used instead of our calculated A^ tot . We can write 
equivalently for our calculation: 



JJ YG f(R,Z,L)dVdL 
^obs / f yD f(R,Z,L)dVdL' 



(16) 
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where f(R,Z,L) is the probability distribution of the pulsar population (§3), Vq and Vo 
are the Galactic volume and the detection volume (for pulsars with pulse period and width 
similar to each of the two observed pulsars), respectively. Note that we do not fix the 
luminosity of the pulsar population to the observed values, but instead, we estimate N tot for 
a distribution of radio luminosities. Since we consider separately PSR B1913+16-like and 
PSR B1534+12-like populations, iV obs = 1. 

For pulsar beaming fractions, we adopt the estimates obtained by KNST: 5.72 for 
PSR B1913+16 and 6.45 for PSR B1534+12. For the lifetime estimates we also follow 
KNST. Our adopted values for the pulsar lifetimes are 3.65 x 10 s yr for PSR B1913+16 and 
2.9 x 10 9 yr for PSR B1534+12. 

Using eq. (15) we calculate -P(72.) for each of the two observed pulsars: 

P(72) = P(N tot " dNtot 



dK 

^WW*. (17) 



5.2. The Total Galactic Coalescence Rate 

Once the probability distributions of the rate contributions of the two observed pulsars 
are calculated, we can obtain the distribution functions of the total coalescence rate 72 to t- 
We define the following two coefficients for each observed system: 

A = (^p) and B = (^*) (18) 

^ Jb ' 1913 ^ Jb ' 1534 

and rewrite the coalescence rate for each binary system: 

P 1913 (K) = A 2 7Ze- An and Pi 534 (^) = B 2 7Ze' Bn . (19) 
One can confirm that each distribution function satisfies the normalization condition 

P(K)dK = 1. (20) 



I 

Jo 



We then define two new variables: 1Z + = 7^1913 + 7^1534 and 1Z- = 7^i 9 i 3 — 7^1534. Since 
72-1913 > and 72i534 > 0, we have — 72 + < 72_ < 72 + . For convenience, we rename 72i 9 i 3 = 
72-1 and 72i534 = 722 and perform a two-dimensional probability distribution transformation: 



P(72 + ,72_) = P(72i,72 2 ) 



dllx 


dK 2 


dn + 


dll- 




dTZ 2 


dK- 


dTZ + 



\P{K U K 2 ). (21) 
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The probability distribution of the total rate TZ to t = Tl+ is obtained after integrating 
P(K + ,Tl-) over 



p(k+)= [ p(n + ,n-) dii- = \ [ p(n 1 ,n 2 )dTZ-. 

Jtl- 2 J n _ 



(22) 



Since the probability distributions for the rate contributions of each of the two observed pul- 
sars are independent of each other, their two-dimensional distribution is simply the product 
of the two. 

p{n l ,n 2 ) = p{n l )P{n 2 ). (23) 

After we replace all variables to 1Z+ and 1Z-, we have 



A 2 B 2 

p(n + = n tot ) = —^e 



2 3 

AB \2 
B - A 



TZ + {e 



J -71+ 



+ e 



B — A 



) • (24) 



We have confirmed that the above function satisfies the normalization / °° P (1Z tot )d1Z tot = 1. 

Having calculated the probability distribution of the Galactic coalescence rate, we can 
take one step further and also calculate ranges of values for the rate 7^ ot at various confidence 
levels (CL). The lower (1Z a ) and upper (7Zb) limits to these ranges are calculated using: 



p(n tot )dn tot = cl 



and 



P(lZ a )=P(1Z h ). 

In all our results, we quote coalescence rates for 68%, 95% and 99% CL. 



(25) 
(26) 



5.3. The Detection Rate for LIGO 

NS-NS binary systems are expected to emit strong gravitational waves during their 
inspiral phase, the late stages of which may be detected by the ground-based gravitational 
wave detectors. In this paper we estimate expected detection rates for LIGO (initial and 
advanced) using the derived probability distribution of the Galactic rate and its values at the 
maximum of the distribution, for different physical models of pulsar properties. To calculate 
the detection rate, we need to extrapolate the Galactic rate to the volume detectable by 
LIGO. We use the ratio between the B-band luminosity density of the Universe and the 
B-band luminosity of our Galaxy as the scaling factor (Phinney 1991; KNST). This is based 
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on the assumptions that (i) the B-band luminosity (corrected for dust absorption) correlates 
with the star-formation rate in the nearby universe and hence the coalescence rate, (ii) the 
B-band luminosity density is constant in the nearby universe. The detection rate, TZdet, is 
calculated by the following equation: 

ftdct = e^tot^det, (27) 

where e is the scaling factor assumed to be ~ 1CT 2 Mpc -3 (for details see KNST). Vdet is 
the detection volume defined as a sphere with a radius equals to the maximum detection 
distance D max for the initial (~ 20Mpc) and advanced LIGO (~ 350 Mpc) (Finn 2001). 



6. RESULTS 

We have calculated the probability distribution of the Galactic coalescence rate 1Z to t, 
its most likely value 7^ pe ak and ranges at different statistical confidence levels, and the most 
likely expected detection rates for the initial and advance LIGO for a large number of model 
pulsar population properties (see Table 2). We have chosen one of them to be our reference 
model based on the statistical analyses and results presented by Cordes & Chernoff (1997). 

For our reference model (L min = 1.0 mJy kpc 2 , p = 2.0, R Q = 4.0 kpc, and Z = 1.5 kpc), 
we find the most likely value of N tot , to be ~ 390 pulsars for the "PSR B1913+16-like" 
population, and ~ 350 pulsars for the "PSR B1534+12-like" population. Using eq. (24), we 
evaluate the total Galactic coalescence rate of NS-NS binaries for this reference case (model 
1 in Table 2). The most likely value of the coalescence rate is ~ SMyr" 1 and the ranges at 
different statistical confidence levels are: ~ 3 — 20 Myr -1 at 68%, ~ 1 — 30 Myr^ 1 at 95%, 
and ~ 0.7 - 40 Myr" 1 at 99%. 

In Figure 4, P(TZ tot ) along with P(7?-i9i3) and ^(7^1534) are shown for the reference 
model. It is evident that the total rate distribution is dominated by that of PSR B1913+16. 
At first this appears to be in contradiction to most other studies of the NS-NS coalescence 
rate (Narayan et al. 1991; Phinney 1991; Curran & Lorimer 1995; van den Heuvel & Lorimer 
1996; Stairs et al. 1998; Arzoumanian et al. 1999; KNST). However, it turns out that this 
difference is due to the fact that earlier studies restricted the scale factor calculation to 
the actual observed pulsar luminosity and the rate estimates were dominated by the low 
luminosity of, and hence large scale factor for, PSR B1534+12. Any corrections to the rate 
estimate which take into account the range of pulsar luminosities were applied as subsequent 
upward corrections. However, this effect is eliminated in our case because we calculate the 
two rate contributions having relaxed the luminosity constraint and instead having allowed 
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for a range in luminosity for the pulsars. In this case any differences in the two separate 
rate contributions depend only on differences in pulse periods, and widths. Given that the 
latter are rather small, it makes sense that, for example, the most likely values of N tot for 
the two pulsars come out to be very similar (e.g. ~ 390 and ~ 350, for PSR B1913+16 
and PSR B1534+12, respectively, in model 1). Consequently any difference in the rate 
contributions from the two populations is due to the difference in lifetimes (about a factor of 
10) for the two observed pulsars (note that the two do not only have similar N tot estimates, 
but also similar beaming correction factors). Since the lifetime estimate for PSR B1913+16 
is much smaller, the total rate distribution is dominated by its contribution. 

In Table 2 we list the population parameters and results of Galactic coalescence and 
LIGO detection rates for a large number of models (but still a subset of the models we have 
investigated in detail). Model 1 is our reference model defined above and the rest of the 
models are used to explore the sensitivity of our results on the pulsar population properties 
in luminosity and space distributions. The variation of the luminosity function parameters 
is within ranges that correspond to 68% confidence levels from the statistical analysis of 
Cordes & Chernoff (1997). The variation of the space scale lengths i?o and Z extends to 
values that are not favored by our current understanding, only to allow us to examine the 
presence of any correlations. Scrutiny of the results presented in Table 2 reveals that rate 
estimates are modestly sensitive (and in some cases, e.g., Z , insensitive) to most of these 
variations, except in cases of very low values of L min , down to 0.3 mJy kpc 2 , and unphysically 
low values of R , down to 2 — 3 kpc. Assumptions about the shape of the space distribution 
(exponential or Gaussian) are not important. The most important model parameter seems 
to be the slope of the luminosity function. The most likely values for the rates are found 
in the range ~ 3 — 22Myr _1 for all models with luminosity-function parameters consistent 
with a 68% confidence level estimated by Cordes & Chernoff (1997) and Rq in the range 
4 — 8 kpc. Within a given model the range of estimated values at 68% confidence level is 
typically broad by a factor of 5, while at the 95% level the ranges broaden by another factor 
of about 5, reaching an uncertainty of ~ 25 for the rate estimates. 

In terms of qualitative variations, it is clear that models with increasing fraction of 
low-luminosity or very distant pulsars lead to increasing rate estimates, as expected. A 
demonstration of this kind of dependence is shown in Figures 5 and 6. We have found that 
there is a strong correlation between the peak value of the total Galactic coalescence rate 
^peak estimated by eq. (24) and the cut-off luminosity L min and the power index p. As seen 
in these Figures TZ pea k increases rapidly with decreasing L min (Fig. 5) or with increasing p 
(Fig. 6). Scale lengths of the spatial distribution (either R or Z ) do not show a correlation 
with 7£ pea k as strong as that of the luminosity-function parameters. In Figure 7 the most 
likely rate estimate is plotted as a function of the radial scale length Rq. It is evident that 
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the rate is relatively insensitive to Ro variations unless Ro becomes very small (< 3kpc). 
Small values of the radial scale for the Galactic distribution imply a large fraction of very 
distant (and hence undetectable) pulsars in the Galaxy and lead to a rapid increase of the 
estimated coalescence rate. However, such small Ro values are not consistent with our current 
understanding of stellar populations in the Galaxy (Lyen et al. 1985). 

7. DISCUSSION 

In this paper we present a new method of estimating the total number of pulsars in 
our Galaxy and we apply it to the calculation of the coalescence rate of double neutron star 
systems in the Galactic field. The method implicitly takes into account the small number of 
pulsars in the observed double-neutron-star sample as well as their distribution in luminosity 
and space in the Galaxy. The modeling of pulsar survey selection effects is formulated in 
a "forward" way, by populating the Galaxy with model pulsar populations and calculating 
the likelihood of the real observed sample. This is in contrast to the "inverse" way of the 
calculation of scale factors used in previous studies. The formulation presented here allows 
us to: (a) calculate the probability distribution of coalescence rates; (b) assign statistical 
significance to these estimates; and (c) quantify the uncertainties associated with them. 

As originally shown by KNST, the most important uncertainties originate from the 
combination of a small-number observed pulsar and a pulsar population dominated by faint 
objects. The probability distribution covers more than 2 orders of magnitude in agreement 
with the uncertainties in excess of two orders of magnitudes asserted by KNST. However, 
for the reference model, even at high statistical confidence level (99%) the uncertainty is 
reduced to a factor of ~ 60. At confidence levels of 95% and 68%, the uncertainty is further 
reduced to just ~ 25 and ~ 5, respectively. We use our results to estimate the expected 
detection rates for ground-based interferometers, such as LIGO. The most likely values are 
found in the range ~ (1 — 30) x 10 _3 yr~ 1 and ~ 4 — 140 yr -1 , for the initial and advanced 
LIGO. 

The statistical method developed here can be further extended to account for distribu- 
tions of pulsar populations in pulse periods, widths, and orbital periods. Most importantly 
the method can be applied to any type of pulsar population with appropriate modifications 
of the modeling of survey selection effects. Currently we are working on assessing the con- 
tribution of double neutron stars formed in globular clusters as well as the formation rate 
of binary pulsars with white dwarf companions that are important for gravitational- wave 
detection by LISA, the space-based interferometer planned by NASA and ESA for the end 
of this decade. 
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Table 1. Simulated Pulsar Surveys 



Telescope 


Year 


v 1 


Av 2 


t 3 

obs 


t 4 

samp 


S 5 ■ 

min 


Detected 6 


Refs 7 


Lovell 76 m 


1972 


408 


4 


660 


40 


10 


51/31 


1,2 


Arecibo 305 m 


1974 


430 


8 


137 


17 


i 


50/40 


3,4 


Molonglo 


1977 


408 


4 


45 


20 


10 


224/155 


5 


Green Bank 300' 


1977 


400 


16 


138 


17 


10 


50/23 


6,7 


Green Bank 300' 


1982 


390 


16 


138 


17 


2 


83/34 


8 


Green Bank 300' 


1983 


390 


8 


132 


2 


5 


87/20 


9 


Lovell 76 m 


1983 


1400 


40 


524 


2 


1 


61/40 


10 


Arecibo 305 m 


1984 


430 


1 


40 


0.3 


3 


24/5 


9 


Molonglo 


1985 


843 


3 


132 


0.5 


8 


10/1 


11 


Arecibo 305 m 


1987 


430 


10 


68 


0.5 


1 


61/24 


12 


Parkes 64 m 


1988 


1520 


320 


150 


0.3 


1 


100/46 


13 


Arecibo 305 m 


1990 


430 


10 


40 


0.5 


2 


2/2 


14 


Parkes 64 m 


1992 


430 


32 


168 


0.3 


3 


298/101 


15,16 


Arecibo 305 m 


1993 


430 


10 


40 


0.5 


1 


56/90 


17-20 


Lovell 76 m 


1994 


411 


8 


315 


0.3 


5 


5/1 


21 


Green Bank 140' 


1995 


370 


40 


134 


0.3 


8 


84/8 


22 


Parkes 64 m 


1998 


1374 


288 


265 


0.1 


0.5 


69/170 


23 


Parkes 64 m 


1998 


1374 


288 


2100 


0.3 


0.2 


-900/600 


24,25 



1 Center frequency in MHz. 
2 Bandwidth in MHz. 
3 Intcgration time in seconds. 
4 Sampling time in milliseconds. 

5 Sensitivity limit in mJy at the survey frequency for long-period pulsars (calculated for each 
trial in the simulations). 

6 Total number of detections and new pulsars 

7 References: (1,2) Davies, Lyne & Seiradakis (1972,3); (3,4) Hulse & Taylor (1974,5); (5) 
Manchester et al. (1978); (6,7) Damasheket al. (1978,1982); (8) Dewey et al. (1985); (9) Stokes 
et al. (1986); (10) Clifton et al. (1992); (11) D'Amico et al. (1988); (12) Nice, Taylor & Fruchtcr 
(1995); (13) Johnston et al. (1992); (14) Wolszczan (1991); (15) Manchester et al. (1996); (16) 
Lyne et al. (1998); (17) Ray et al. (1996); (18) Camilo et al. (1996); (19) Foster et al. (1995); 
(20) Lundgren, Zcpka & Cordes (1995); (21) Nicastro et al. (1995); (22) Sayer, Nice & Taylor 
(1997); (23) Edwards et al. (2001); (24) Lyne et al. (2000); (25) Manchester et al. (2002) 
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Table 2. Model Parameters and Estimates for lZ to t and IZdet & t Various Statistical 
Confidence Levels for Different Population Models 



Model 1 


Parameters 




fttot (Myr- 1 ) K det of LIGO (yr" 1 ) 


T2 
min 


p 3 Ri 




peak 7 68% 8 95% 8 initial (xl0~ 3 ) advanced 


(mJy kpc 2 ) 


(kpc) 


(kpc) 


peak 7 68% 8 peak 7 68% 8 



1 


1.0 


2.0 


4.0 


(G 6 ) 


1.5 


(E 6 ) 


8.0 


2 


1.0 


2.0 


4.0 


(G) 


0.5 


(E) 


7.1 


3 


1 n 

1.0 


z.O 


4.0 


(G) 


z.O 


(E) 


O A 

8.4 


4 


1.0 


Z.O 


4.0 


( E ) 


1.5 


(E) 


8.7 


r- 




1.0 


Z.O 


A n 

4.0 


(G) 


1.5 


(G) 


7.9 


6 


0.3 


2.0 


4.0 


(G) 


1.5 


(E) 


26.9 


7 


0.7 


2.0 


4.0 


(G) 


1.5 


(E) 


11.5 


8 


1.5 


2.0 


4.0 


(G) 


1.5 


(E) 


5.5 


9 


3.0 


2.0 


4.0 


(G) 


1.5 


(E) 


2.9 


i n 
1U 


n q 
U.o 


l.O 


a n 
4.U 




1.0 




O A 

y.4 


11 


0.7 


1 o 
1.0 


A n 

4.0 


(G) 


1.5 


(E) 


A O 

4.8 


12 


1.0 


1 o 

1.8 


A n 

4.0 


(G) 


1.5 


(E) 


3.6 


lo 


1.0 


l.O 


a n 
4.U 


(G) 


1.0 


(E) 


z. 1 


1 /I 

14 


o.U 


1.0 


a n 
4.U 


(Gj 


1.0 


(EJ 


1.0 


15 


0.3 


2.2 


4.0 


(G) 


1.5 


(E) 


61.2 


16 


0.7 


2.2 


4.0 


(G) 


1.5 


(E) 


22.1 


17 


1.0 


2.2 


4.0 


(G) 


1.5 


(E) 


14.9 


18 


1.5 


2.2 


4.0 


(G) 


1.5 


(E) 


9.8 


19 


3.0 


2.2 


4.0 


(G) 


1.5 


(E) 


4.7 


20 


1.0 


2.5 


4.0 


(G) 


1.5 


(E) 


28.3 


21 


1.0 


2.0 


2.0 


(G) 


1.5 


(E) 


26.1 


22 


1.0 


2.0 


3.0 


(G) 


1.5 


(E) 


12.8 


23 


1.0 


2.0 


5.0 


(G) 


1.5 


(E) 


6.7 


24 


1.0 


2.0 


6.0 


(G) 


1.5 


(E) 


6.6 


25 


1.0 


2.0 


7.0 


(G) 


1.5 


(E) 


6.9 


26 


1.0 


2.0 


8.0 


(G) 


1.5 


(E) 


7.4 



+9.3 

-4.7 

+8.3 

-4.2 

+9.9 

-5.0 

+10.2 

-5.1 

+9.2 

-4.6 



+32.0 

-16.1 

+13.5 

-6.8 

+6.3 

-3.2 

+3.3 

-1.7 



+10.8 

-5.5 

+5.4 

-2.8 

+4.1 

-2.1 

+3.0 

-1.5 

+1.8 

-0.9 



+75.8 

-37.5 

+27.0 

-13.4 

+18.0 

-9.0 

+11.7 

-5.9 

+5.5 

-2.8 



+35.6 
-17.5 



+29.7 

-15.2 

+14.6 

-7.4 

+7.9 

-4.0 

+7.8 

-3.9 

+8.2 

-4.1 

+8.8 

-4.4 



+23.3 

-6.7 

+20.8 

-5.9 

+24.7 

-7.1 

+25.6 

-7.3 

+23.0 

-6.6 



+80.3 

-22.7 

+33.9 

-9.7 

+15.9 

-4.6 

+8.3 

-2.4 



+27.1 

-7.9 

+13.5 

-4.0 

+10.3 

-3.0 

+7.6 

-2.2 

+4.4 

-1.3 



+190.3 

-52.1 

+67.8 

-18.7 

+45.2 

-12.6 

+29.4 

-8.2 

+13.8 

-3.9 



+89.4 
-24.2 



+74.2 

-21.7 

+36.4 

-10.7 

+19.8 

-5.6 

+19.5 

-5.5 

+20.5 

-5.8 

+22.2 

-6.3 



3.3 
3.0 
3.5 
3.6 
3.3 

11.3 
4.8 
2.3 
1.2 

3.9 
2.0 
1.5 
1.1 
0.7 

25.6 
9.2 
6.2 
4.1 
2.0 

11.8 

10.9 
5.4 
2.8 
2.7 
2.9 
3.1 



+3.9 
-2.0 
+3.5 
-1.7 
+4.1 
-2.1 
+4.3 
-2.2 
+3.9 
-1.9 



+13.4 

-6.7 

+5.7 

-2.9 

+2.7 

-1.3 

+1.4 

-0.7 



+4.5 
-2.3 
+2.3 
-1.2 
+1.7 
-0.9 
+1.3 
-0.6 
+0.7 
-0.4 



+31.7 

-15.7 

+11.3 

-5.6 

+7.5 

-3.8 

+4.9 

-2.5 

+2.3 

-1.2 



+14.9 
-7.3 



+12.4 

-6.3 

+6.1 

-3.1 

+3.3 

-1.7 

+3.3 

-1.6 

+3.4 

-1.7 

+3.7 

-1.9 



17.9 
15.9 
19.0 
19.5 
17.7 

60.5 
25.9 
12.3 
6.4 

21.2 
10.7 
8.2 
6.0 
3.5 

137.6 
49.7 
33.5 
22.0 
10.5 

63.6 

58.6 
28.9 
15.1 
14.8 
15.5 
16.8 



+21.0 

-10.6 

+18.7 

-9.4 

+22.2 

-11.2 

+23.0 

-11.6 

+20.7 

-10. 



+72.1 

-36.2 

+30.5 

-15.3 

+14.3 

-7.2 

+7.4 

-3.8 



+24.3 

-12.4 

+12.2 

-6.2 

+9.3 

-4.7 

+6.8 

-3.5 

+4.0 

-2.0 



+170.5 

-84.4 

+60.8 

-30.2 

+40.5 

-20.2 

+26.4 

-13.2 

+12.4 

-6.2 



+80.0 
-39.4 



+66.7 

-34.1 

+32.8 

-16.8 

+17.8 

-8.9 

+17.5 

-8.8 

+18.4 

-9.2 

+19.9 

-10.0 
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Table 2 — Continued 



Model 1 




Parameters 




n tot (Myr- 


- 1 ) 


K dct of LIGO (yr- 1 ) 




7-2 
min 


n 3 R 4 
p 2t 




peak 7 68% 8 


95% 8 


initial (xlO -3 ) advanced 




(mJy kpc 2 ) 


(kpc) 


(kpc) 






peak 7 68% 8 peak 7 68% 8 


27 


1.0 


2.0 9.0 (G) 


1.5 (E) 


8 - 4 -5.0 


+25.1 
-7.1 


3.5 t\\ 18.9 ±g : | 



iModel No. 

2 Minimum luminosity L m i n in mJy kpc 2 . 
3 Power index of the luminosity function p. 
4 Radial scale length R in kpc. 
5 Vertical scale height Z G in kpc. 

6 Gaussian (G), and exponential (E) functions for spatial distributions. 
7 Peak value from P(R to t)- 
8 Confidence level. 



-22 - 




Integration time (min) 



Fig. 1. — Average signal-to-noise degradation factor in pulsar search code versus survey 
integration time for PSR B1913+16 and PSR B1534+12. 
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Fig. 2. — The Poisson-distribution fits of P(N obs ) for three values of the total number N tot 
of PSR B1913+16-like pulsars in the Galaxy (results shown for model 1). Points and error 
bars represent the counts of model samples in our calculation. Dotted lines represent the 
theoretical Poisson distribution. 
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Fig. 3. — The linear correlation between A =< iV obs > and N tot is shown for model 1. 
Solid and dashed lines are best-fit lines for PSR B1913+16-like and PSR B1534+12-like 
populations, respectively. Points and error bars represent the best-fit values of A for different 
values of A^tot- 
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Fig. 4. — The probability distribution function of coalescence rates in both a logarithmic and 
a linear scale (small panel) is shown for model 1. The solid line represents P(1Z tot ) and the 
long and short dashed lines represent P(H) for PSR B1913+16-like and PSR B1534+12-like 
populations, respectively. We also indicate the confidence levels for P(1Z tot ) by dotted lines. 
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Fig. 5. — The correlation between 7£ pca k and the cut-off luminosity L min with different power 
indices p of the luminosity distribution function. 
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Fig. 6. — The correlation between 7£ pca k and the power index of the luminosity distribution 
function p. 
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Fig. 7. — The correlation between 7£ pca k and the radial scale length R . 7^ pea k is not sensitive 
to R in the range between 4-9 kpc. 



